Collapse of a Bose-Einstein condensate induced by fluctuations of the laser intensity 
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M-H The dynamics of a metastable attractive Bose-Einstein condensate trapped by a system of laser 

beams is analyzed in the presence of small fluctuations of the laser intensity. It is shown that 
the condensate will eventually collapse. The expected collapse time is inversely proportional to 
the integrated covariance of the time autocorrelation function of the laser intensity and it decays 
logarithmically with the number of atoms. Numerical simulations of the stochastic 3D Gross- 
Pitaevskii equation confirms analytical predictions for small and moderate values of mean field 
interaction. 



PACS numbers: 03.75.Kk, 42.65.-k, 42.50.Ar 



I. INTRODUCTION 



. The experimental realization of Bose-Einstein condensation (BEC) in dilute atomic gases founded a rapidly 

fSJ ' progressing new field of research Q. The physical properties of BECs, which to date comprise eight elements Rb, Na, 
£NJ , Li, H, He, K, Cs, Yb and their isotopes, are predominantly determined by interatomic forces. Some of the atomic 
species ( 7 Li, 85 Rb, 133 Cs) possess a negative s-wave scattering length in the ground state and display attractive 
interactions. The attractive interaction between the atoms causes the collapse of the BEC so that a stable BEC was 
not believed to exist H- However, when an external spatial confinement is imposed for instance by a system of laser 
beams, a trapping potential shows up which can counterbalance the attractive interaction and allows the formation of 
a metastable BEC. When the number of atoms increases, the attractive interaction becomes stronger and the energy 
barrier that prevents the 3D BEC from collapsing becomes weaker. To a given trapping potential there corresponds a 
critical number of atoms above which the energy barrier vanishes. The case of a quadratic potential has been studied, 
the critical number of atoms has been computed by a variational approach and by extensive numerical simulations of 
Sh ! the Gross-Pitaevskii (GP) equation, and the results have been checked experimentally 0,0,0- 

One of the most important aspects of BECs in the regime of attractive interactions is that they are unstable against 
collapse. The collapse shows up as a rapid and strong shrinking of the condensate at some critical number of atoms, 
and is accompanied by significant atomic losses due to many-body processes 8j. The collapse is initiated when the 
balance of forces governing the size and shape of the condensate is altered either by internal or external factors. With 
respect to spatial and energetic stability the magnetic traps appear to be better controllable compared to optical 
traps On the other hand, due to increasing interest in far-off resonant laser traps for Bose-condensation of atoms 
which are insensitive to magnetic fields [lfj , the investigation of different aspects of BEC dynamics in optical traps 
is becoming a very relevant subject. Of particular interest is the effect of temporal fluctuations of the laser intensity 
which in turn involve temporal fluctuations of the parabolic trapping potential |hJ . In the present paper we shall 
consider the BEC dynamics under random fluctuations of the strength of the parabolic trap potential and we shall 
show that small fluctuations can lead to the eventual collapse of the 3D BEC due to a cumulative effect of stochastic 
perturbations. The random fluctuations have all harmonics in their spectrum, and some of them participate in the 
parametric resonance leading to collapse. This stochastic parametric resonance in the BEC width oscillations has a 
rough equivalent particle picture: the Kramers' exit problem which is concerned with noise activated escape from a 
potential well |l2| . 

Quantum tunneling (QT) is considered as playing a key role in the condensate collapse when the number of atoms is 
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close to the critical number [13j . We shall see that the BEC instability driven by random fluctuations of the strength 
of the parabolic trap potential is all the more dramatic as the number of atoms is closer to the critical number. 
Our consideration thus shows that even weak noise can play a competitive role in this limit with QT and should be 
taken into account. The effect of optical trap noise was previously considered in the context of stochastic heating 
of trapped atoms [Til IT^ . In a far-off resonant optical trap created by a system of red detuned lasers the variable 
trapping potential can be represented as V(t,r) = — a\E(t, r)| 2 /4, where a is the atomic polarizability and E(t,r) 
is the electric field amplitude. The dynamics of trapped atoms can be described by the corresponding Hamiltonian 
H = p 2 /(2m) + (1/2)to(jJq(1 + T](t))r 2 , where uiq = k^/m is the mean square trap oscillation frequency, and fco is 
proportional to the time- averaged laser intensity Iq ~ \E\ 2 . The time dependent spring constant is determined by 
fractional fluctuations of the laser intensity rj(t) = (I(t) — Iq)/Iq ^3 ■ The influence of the fluctuations of the 
trap potential on the dynamics of ID GP type equation has been considered in [l5| and the trap and nonlinearity 
fluctuations in two dimensional BEC in 

The paper is organized as follows. In Section ^ we give a description of the model and apply a variational 
approach. In Section ITTT1 we derive the effective dynamics of the action-angle variables of the system driven by random 
perturbations. Section lTvl frcsp. Ivjl are devoted to the asymptotic analysis of the system for small (resp. near-critical) 
number of atoms. Finally we check the variational approach and our asymptotic analysis in Section IVTI bv performing 
direct numerical simulations of the GP equation. 



II. THE MODEL AND THE VARIATIONAL APPROACH 

We consider the mean-field GP equation for the single-particle wave function 0] 

h 2 

th4> t = - — AiP + V{t,r)iP + g\^\ 2 ^. (1) 
2m 

The nonlinear coefficient is g = 4irh 2 a s /m where a s and m are respectively the atomic scattering length and mass. 
The number of atoms is N = J \i[>\ dx. V is the external trapping potential imposed by a system of laser beams. 
We consider a harmonic model, but we take into account temporal fluctuations of the laser intensity which in turn 
induces temporal fluctuations of the quadratic potential 



muj, 



V(t,r) = ^\v\ 2 [l + r,(t)]. (2) 

For the optical trap oj 2 — al / (2ml 2 ,) , where Iq is the size of the laser beam, / is the intensity, a is a constant 
proportional to the laser frequency detuning. The random function r)(t) describes the laser intensity fluctuations 
r](t) = (I(t) — Iq)/ Iq. The stationary random process ry has zero-mean and standard deviation a n . We shall see in the 
following that the standard deviation is not sufficient to predict the collapse of the BEC, but the coherence time and 
more generally the power spectral density of rj will play a role. 

We now cast Eq. Q in a dimensionless form by introducing the variables t' — mot, r' = r/ro, r^ 1 — \J muio/h, and 
u = y/4Tr\a s \rQip. This yields the following partial differential equation (PDE) 

iu v = -\b!u + i|rf [1 + r,'(t')]u + a s \u\ 2 u, (3) 

where a s = sgn(a s ) = ±1 and rj'(t') — rj(t' /wq). From now on we drop the primes. The next step consists in applying 
the variational approach. This approximation was first introduced by Anderson |l9j | and developed in nonlinear optics 
pof . A similar technique was elaborated for the BEC dynamics based on the GP equation [2l| . The variational ansatz 
for the wave function of the BEC is chosen as the Gaussian Q 

u(t, r) = A(t) exp (-^L + l M?L + «, (f) ) . (4 ) 



a 



(0)ro is the initial BEC rms width in physical variables 



V2 

The number of atoms is 



( / |r| 2 |^ = 0,r)| 2 d 3 r 



1/2 
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FIG. 1: Potential U(a) for P — 0.2. The important points (ai < ao < 02) are also represented. 



Following the standard procedure j2£|, we substitute the ansatz into the Lagrangian density generating Eq.© and 
calculate the effective Lagrangian density in terms of A, a, 6, 8 and their time-derivatives. The evolution equations 
for the parameters of the ansatz are then derived from the effective Lagrangian by using the corresponding Eulcr- 
Lagrange equations. In particular this approach yields a closed-form ordinary differential equation (ODE) for the 
BEC width a 

Ott + a(l + i 7 (t)) = i + ^V ) (5) 

where P = yj2/irN\a s \/rQ. We study in this paper the attractive case (a s < 0, <x s = —1). The evolution equation 
finally reads 

a u +a(l + r](t)) = ^-^. (6) 



III. ACTION-ANGLE VARIABLES 
A. Unperturbed dynamics 

The energy E of the unperturbed BEC is given by: 

E(t) = ±a 2 t (t) + U(a(t)), U {a)= 1 -(a 2 + ^j-^. (7) 

In absence of random fluctuations 77 = the energy E is an integral of motion. The BEC width obeys a simple 
dynamics with Hamiltonian structure 

H(p,q) = ±p 2 + U(q) (8) 

with q = a and p = a t . A straightforward analysis 0,H3 shows that if P < P c = 4/5 5 / 4 ~ 0.535, then the potential 
U possesses a local minimum that we shall denote by ao (see Fig. QJ. The corresponding ground state has energy 
Eq = [/(ao). Below ao there is the local maximum ai with energy E\ = U(a±), and below a\ the potential decays to 
—00. Above ao the potential increases to +00. It crosses the energy level E± at 122 • 

If the initial conditions (a(0),a t (0)) correspond to an energy above Ex, or below Ex but a(0) < ax, then the 
condensate width goes to zero in finite time which means that the BEC collapses. On the contrary, if the initial 
conditions (a(0), a t (0)) correspond to an energy between Eq and Ex, and a(0) > ai, then the orbits of the motion are 
closed, corresponding to periodic oscillations. In order to explicit the periodic structure of the variables a and a t , we 
introduce the action-angle variables. The orbits are determined by the energy imposed by the initial conditions: 



E=-a 2 t (Q) + U(a(0)). 
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For E G (Eo, Ei), we introduce 6\(E) < 62(E) the extremities of the orbit of a for the energy E: 

U(e 1 (E)) = U(e 2 (E))=E. 
The action / is defined as a function of the energy E by 

r e 2 {E) 



1 r 1 pe^yaj 

1(E) = — <ppdq= - I ^2E - 2U(b)db. 

27T / 7T J ei (E) 



(9) 



The motion described by JBj is periodic, with period 



T(E)=^ = < i r m / db (10) 



or else T(£') = 2ir^(E). The angle is defined as a function of / and a by 

r <9P , 2vr f a db 

<j)(E, a) = - -£dq ' 



81 T(E) J y/2E-2U(b)' 

The transformation (E, a) — > (J, 0) can be inverted to give the functions £ (J) and -4(J, 0). The BEC width oscillates 
between the minimum value e\(E) and the maximum value 62(E). The energy E as well as the action / are constant 
and fixed by the initial conditions, so the evolution of the BEC width is governed by 

m = m - ^Jy*. 

B. Perturbed dynamics 

From now on we assume r\ ^ and we denote by <J n the standard deviation of r\. We investigate the stability of 
the BEC when the unperturbed motion is oscillatory. In particular we aim at studying the collapse time T c defined 
as the first time t such that a(t) = 0. While the energy of the BEC is below Ei, the orbit is closed. As soon as the 
energy reaches the energy level Ei, the BEC collapses in a time of order 1 (w.r.t. a„). We shall show that the hitting 
time for the energy level E\ is of order cr~ 2 , so the collapse time T c is imposed by the hitting time Th defined as the 
first time t such that E(t) = E\ or equivalently I(t) = I\ :— T(Ei). 

In presence of perturbations, the motion of a is not purely oscillatory, because the energy and the action are slowly 
varying in time. We adopt the action-angle formalism, because it allows us to separate the fast scale of the locally 
periodic motion and the slow scale of the evolution of the action. Thus, after rescaling r = a^t the action-angle 
variables satisfy the differential equations 

°n u v Mil 
— = 2 W ( J ) ri( — )hi(I,<P)i 

where h(I,(f)) = ^A 2 (I,4>) and oj(I) — YffiJS) are sn^ooth functions and h is periodic with respect to cj) with period 
2ir. The normalization r = a ^t h as been chosen so that the random process 77 appears with the scales of a white noise 
in the differential equations (|llf> . Applying a standard diffusion-approximation theorem |23j . we get that (I(t)) t >o 
behaves like a diffusion Markov process with the infinitesimal generator 

1 d 2 d 

where 

,-2tt 



-1 />Z7T /"OO 

A(I) = - / h^I, <j>)h4l, <\> ~ uj(I)t)E[ri(0)ri(t)}dtdct>, 

n Jo Jo 

-1 p2lT f>00 

B(I) = - / ht(I,<l>)hv(I,<l>-Lj(I)t)E[Ti(0)Ti(t)]dtd(l>. 
k Jo Jo 
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This means in particular that the probability density function of I(t) satisfies the Fokker-Planck equation dtp = £}p, 
pit = } I) = 5(1 — To); where I is the initial action at time and C*j is the adjoint operator of £j, i.e. £*jp = 
(l/2)dj [A(I)p] — di [B(I)p\. Moreover, standard results of stochastic analysis allow us to compute recursively the 
moments of Th |3. Denoting I\ — T{E\), the first moment fjp-^I) — E/[T^] (the mean value of T), starting from 
action I at time 0) satisfies 

£jM (1) = -1, /i (1) (A)=0. (12) 

The rc-th moment ^(1) = E/[T£] satisfies 

= -n^ n - x \ fi (n) (h) = 0. (13) 

In the following sections we shall apply and discuss these general results in two different frameworks: small and critical 
nonlinearity. 



IV. SMALL NONLINEARITY 



A. Expansions of the action-angle variables for small nonlinearity 



In this section we assume that P -C 1 which will allow us to derive simple expressions for the physically relevant 
quantities. The points a,j and Ej can be expanded for small nonlinearity P as 

a = l + O(P), ai = P + 0{P 2 ), a 2 = J=-+0(l), 

E = l + O(P), £ 1 = ^L + 0(I). 

Note that, as P becomes small, the potential barrier grows like P~ 2 , which shows that the trap looks like a deep 
quadratic external potential. The functions h(I,(f>) and uj(I) can also be expanded for any <f> and I < I\ = I(E\) = 
1/(12P 2 ) + 0(1/P): 

h(I, 4>) = -+I+ V/ + I 2 cos(0) + O(P), 
w(I) = 2 + 0(P). 

Accordingly the unperturbed dynamics of the BEC width for small P is approximately 

a(t) = \jl + 2I + 2^I +I 2 cos(2f). (14) 

Figure [3 shows that this approximation is indeed very good for the orbit a(t) whatever the initial conditions lying in 
a closed orbit with energy < Ex. 



B. Effective equations in presence of perturbations 



In case of small nonlinearity P -C 1, the above expansions allow us to derive simple effective equations for the BEC 
action in presence of perturbations. Applying the general results obtained in Section flll Bl we get that the action I(t) 
behaves like a diffusion process with the infinitesimal generator 



r 1 9 

Cl = 2 ac dl 



(I + I 2 ) 



d_ 

01 



where 



cos(2i)E[i7(0)»7(t)]dt. 
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FIG. 2: Unperturbed dynamics of the BEC width. We assume P = 0.1, a t (0) = 0, a(0) = 2 (picture a), a(0) = 5.7 (picture b). 
The second case corresponds to an energy very close to Ei. The results from the resolution of the ODE are compared with the 
asymptotic formula 1141 . 



The expression of Ci holds true only for I < I±. We can compute the growths of the first moments of the action 
starting from the ground state I — while e Qct <C P~ 2 : 



E Q [/(*)] : 

Mm 2 ] 



a t 



1 

2' 



= 3a c t 



(15) 
(16) 



An empirical way to estimate the mean disintegration time is to look for the time t± such that Eo[/(ii)] = Ji, where 
Ii = 1/(12P 2 ). From Eq. H15[l we get t\ = (l/a c )ln[l + 1/(6P 2 )]. This argument is rough because the expectations 
are ill-placed. The exact results provided by the rigorous stochastic analysis confirm that this prediction is not correct. 
Integrating Eqs. (|12I13|) we get that the expectation of the disintegration time starting from the ground state / = is 



Eo\T h ] = — ln(l + =■) 



p«i 2 



(-21n(P) - ln(12)) , 



(17) 



while its variance is 



Varo(Th) = -j 



F<1 



ln(l 



12P 2 



) + dilog(l + 



1 x 1 , / 

+ - n(l 

12P 27 2 V 



-21n(P)-ln(12)- — 
6 



where the dilogarithm function is the tabulated function defined as follows: 

m (y) 



dilog(x) 



i-y 



dy. 



12P 2 ' 



(18) 



Equations H17I18|I are the most important results of this paper. They show that the collapse time varies as ~ ln(P 2 ), 
while the energy barrier is ~ P -2 . In physical variables, the expected collapse time is 



MTc] = hi 1 

UJqO. 



Tin 



24mLU a 2 s N 2 



Taking the experimental data ujq = 10kHz, N ~ 5 ■ 10 3 , a s 
collapse time « (1 — 10) seconds. 



a = uj cos(2u) t)E[r)(0)r](i)]dt. 
Jo 

- — 5nm, and a = 10 -4 — 10 -5 , we obtain the expected 
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TABLE I: Comparisons between the averages and rms of the collapse time obtained from numerical simulations and from 
theoretical formulas. Here a — 0.3 and t c = 0.5. 



p 


(r) 


rms(r) 


num 


theor 


error 


num 


theor 


error 


0.05 


4112 


4103 


0.2% 


2241 


2335 


4% 


0.1 


2585 


2591 


0.2% 


1718 


1601 


7% 


0.2 


1257 


1306 


3.5% 


833 


865 


4% 


0.3 


586 


760 


23% 


407 


518 


21% 


0.4 


205 


486 


58% 


165 


336 


51% 




FIG. 3: Histograms of the collapse time obtained from series of 1000 simulations. Picture a: P = 0.1. Picture b: P = 0.05. 



C. Numerical simulations 



We compare the theoretical predictions with numerical simulations of the ODE |JB}. We use a fourth-order Runge- 
Kutta method for the resolution of the ODE. The random fluctuations are modeled by a stepwise constant random 
process: 

3 

where the Xj are independent and identically distributed random variables with uniform distribution over (—1/2, 1/2) 
and t c is the coherence time of the laser. The coefficient a c is then given by 

2 1 — cos(2i c ) 



The first series of simulations were performed with the parameters a — 0.3 and t c = 0.5. We investigate different 
configurations corresponding to different values of the parameter P starting from a(0) = 1, a t (0) = which is very 
close to the ground state. We have carried 1000 simulations for each configuration. The theoretical values for the 
expected value and standard deviation according to formulas <|17I18H are reported in Tableland compared with the 
values obtained from averaging of the results of the numerical simulations. 

Note that the statistical formulas are theoretically valid in the asymptotic framework P <C 1. The numerical 
simulations show that they are actually valid for P < 0.2. More exactly, the comparisons between the theoretical 
predictions and the numerical simulations shows excellent agreement for the mean values, and very good agreement 
also for the standard deviations. We also plot in Fig.[3]thc histograms of the collapse times for two series of simulations. 

Finally, in Table ITT1 we report results with a high level of fluctuations (namely a = 2). The theoretical predictions 
are still in agreement with the numerical simulations for P < 0.3 with an accuracy of 10% although the considered 
configurations are at the boundary of the validity of the asymptotic theory. 
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TABLE II: Comparisons between the averages and rms of the collapse time obtained from numerical simulations and from 
theoretical formulas. Here a — 2 and t c = 0.5. 



p 




rms(r) 


num 


theor 


error 


num 


theor 


error 


0.05 


98.6 


92.3 


6.4% 


55.5 


52.5 


5.4% 


0.1 


63.7 


58 


8.5% 


39.1 


36.0 


7.9% 


0.2 


31.9 


29.4 


7.8% 


21.2 


19.5 


8.2% 


0.3 


16.1 


17.1 


6.5% 


11.4 


11.7 


2.6% 


0.4 


6.6 


10.9 


65% 


4.9 


7.6 


55% 



0.762 



(a) 




FIG. 4: Picture a: Potential U(a) for P = P c - 0.01 ~ 0.525 (<$ = 0.01). Picture b: Unperturbed dynamics of the BEC width. 
We assume a t (0) = 0, a(0) = 0.66, 5 = 0.01. The results from the resolution of the ODE are compared with the asymptotic 
formula 



V. CRITICAL NONLINEARITY 



A. Expansions of the action-angle variables for critical nonlinearity 

In this section we address the case where the nonlinear parameter P is close to the critical value P c = 4/5 5 / 4 . We 
do so by setting P = P c — 8 and assuming 5 <C 1. Once again, all quantities can be expanded in powers of 8. After 
some algebra, we get 

aj = a g + 2- 1/2 5~ 1/8 a j( 5 1/2 + 0(6) with a g = 5~ 1/4 , 5 = 1, &i = -1, a 2 = 2, 

Ej =E g + 2 1 /2 3 -i 5 V8£:,(53/2 + 0{8 2 ) with E g = Z-W 2 + 3- l 5 3 '% E Q = -1, ^ - 1. 

More generally, if a € [ai,a 2 ], then it can be parameterized as a = a g + 2~ 1 / 2 5~ 1 / 8 (5 1 / 2 a and the potential at a can 
be expanded as 

U(a) =E g + 2 l 'H- 1 b 7 'H^ 2 U{a) + 0(8 2 ), 

where 

U(~a) = ^(~a 3 -3~a). 

Note that locally (i.e. around a g ) the potential presents a local minimum at ao (see Figure but the shape of the 
potential well is very different from the one observed in the framework P < 1 (compare with Fig. [IJ. The width of 
the well et2 — a\ is of the order y/d and its depth E\ — Eq is of order <5 3 / 2 . The local shape of the potential is given 
by the cubic function U. 

We now consider the action-angle variables. If E € [E\,E2), then it can be parameterized as E = E g + 
2 1 / 2 3- 1 5 7/8 E8 3/2 with E G [E ,Ex). There exist three solutions e 3 (E) < 5i < e^E) < e 2 (E) < 5 2 of the cu- 
bic equation U(a) = E. e\(E) and &2{E) determine the extremities of the orbit of the normalized width a for the 
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normalized energy E in case of unperturbed dynamics. The cubic equation can be solved: 

e i (^) = 2co S ^ CCOS ^+ 2 ^- 2 )^. 

In particular, if E = Eq (i.e. E = Eq), then e\(Eo) = ^{Eq) = 1 (and e3(_Eo) = —2), which corresponds to the 
ground state a(t) = Oo, or a(t) = 1. 

The period T(E) of the closed orbit at energy level E, as defined by (|10|) . can be expanded as well. Introducing 

f(E) = 5 1 /4 2 -V4 5 9/16 T( ^ g + 2 1 /2 3 -l 5 7/8^ (5 3/2 )] 

we get that T is at leading order with respect to S a 0(l)-function that can be expressed in terms of tabulated 
functions 

t{E) = 275 JT (p(£)) , 

^Je 2 (E)~e 3 (E) 

where 

e 2 (S)-ei(S) 



e 2 (£)-e 3 (£) 



and if is the complete elliptic integral [2^, p. 590]. We then define a normalized action T{E) for G [i?o, Ex] = [— 1, 1] 
I'.v 



1 ^ 



- — / r( fl )dfl. 



The function J : [E ,Ei] — > [0,h] is invertible. Its inverse is denoted by £ : [0, 1{\ — > [i?o,-Ei] where Ji = 18/(57r). It 
is plotted in Fig. We can see that £ is roughly linear. Similarly we can define the angle (f>(E, a) and its inverse 
A(I, 4>). The function A : [0, I\] x [0, 2tt) — > [01,03] can be expressed in terms of Jacobian elliptic functions 



M!,*) = ei (£ (I)) + [e 2 (£(I)) - ei (£ (/))] sn 2 f ^ ( ^ (J)) ^p(£(J))j 



(19) 



where sn is the Jacobian sinus [2^ p. 589]. In absence of perturbation the action is preserved and the closed orbit of 
a(t) for a normalized action I £ [0, Ji) is given by 

3(t) = i (/, *(t)) With ^(t) = _^/4 2 -l/4 3 -l/ 25 9/16^^_ (20) 

The true orbit is a(i) = a g + 2~ 1 / 2 5~ 1 / 8 (5 1 / 2 a(t). Figure^J) shows that this approximation (derived in the asymptotic 
framework (5 1) is indeed reasonably good. 

B. Effective equations in presence of perturbations 



Following the strategy presented in Section llll Bl we introduce the normalized action-angle variables so that E(t) = 
£(I(t)) and a(t) = A(I(t),4>(t)). While the energy of the BEC is below Ei, the orbit is closed. As soon as the energy 
reaches the energy level E\, the BEC collapses in a time of order 1 (w.r.t. a v ). We shall show that the hitting time 
for the energy level E\ is of order <r~ 2 , so the collapse time T c is imposed by the hitting time Th defined as the first 

time t such that I(t) = I\. Here we rescale r = cx 2 (5 _3 / 2 i. This normalization is chosen so that the random process r\ 
appears with the scales of a white noise in the differential equations 

(xt e 5 

# = _^V W(/) _i, ( ^ )h{l4>) , 
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(b) 




FIG. 5: Functions 7 i— » £{I) (picture a) and 1 1— ► ^4(7) (picture b). 



where e = a^S^ 4 , h(I,<f>) = 2- 5 / 4 3 5- n / lf U(7», and cD(I) 



Note once again that h and uj are smooth 



T(£(z)) ■ 

functions, and ft. is periodic with respect to </> with period 27r. By applying a diffusion approximation theorem |23j| . 
we get that (I(t))t>a behaves like a diffusion Markov process with the infinitesimal generator 

1 dl dl 



where 



A(I) = 2- 1 / 2 3 2 5- 11 / 8 [e 2 (f(7))- ei {£{!)) 
E[r)(0)r)(t)}dt, 



2 K(p(£(I)) f K (P^ 



cn 2 d 



n 2 sn 2 ^s, p(f(7))J ds, 



and dn and cn are two tabulated elliptic functions |25l p. 589]. The conditions ensuring the diffusion-approximation 
are S <C 1, cr 2 <§; S 3 ^ 2 . The diffusion coefhcient A(I) is plotted in Fig. 03. 

Using the results reported in Section IlII Bl we get the following recursive relation (n > 1) for the moments of the 
hitting time 



11 hi a J i A(x) 



(21) 



where 7i = 18/(57r). In dimensional variables, the result reads as follows. Starting from the ground state ao, the 
expected value of the collapse time is 



(P c -P) 3/2 
^opc = o Cl, 



(22) 



where C\ is the constant Ci = J Q 1 ^^y^- By a numerical integration using Matlab we have found Ci ~ 8.5. More 
generally, we have 

(P c - P) 3n/2 



u 2n a 



~Cni 



(23) 



where C n are constants obtained recursively from Eq. I|21|) . By a numerical integration we have found Ci ~ 110. 



C. Numerical simulations 



We compare the theoretical predictions with numerical simulations of the ODE JJJJ. We use the same model as 
in Section TlV CI with the parameters a = 0.025 and t c — 0.5. We report in Table ITU the theoretical values for the 
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TABLE III: Comparisons between the averages and rms of the collapse time obtained from numerical simulations and from 
theoretical formulas. 



p 


8 


(r) 


rms(r) 


num 


theor 


error 


num 


theor 


error 


0.525 


0.01 


651 


653 


0.3% 


447 


472 


5% 


0.515 


0.02 


1754 


1846 


5% 


1240 


1334 


7.5% 


0.505 


0.03 


3175 


3392 


7% 


2217 


2451 


10.5% 


0.495 


0.04 


4673 


5222 


11.5% 


3107 


3775 


21.5% 



expected value and standard deviation according to formulas <|22I23[1 as well as the values obtained from averaging of 
the results of the numerical simulations. The statistical formulas are theoretically valid in the asymptotic framework 
S(= P c — P) <C 1. The numerical simulations show that they are actually valid for S < 0.03. 

VI. VALIDATION OF THE VARIATIONAL APPROACH 

The analysis carried out in this paper is based on the variational approach using a Gaussian ansatz. The Gaussian 
ansatz for the study of static and dynamic properties of trapped gases has been widely used (see for instance |2ll l26l 
113, The variational principle is shown in these papers to be a simple Lagrangian-based method that gives 

reasonable accurate ordinary differential equations approximations to the true dynamics for the solution of the GP 
equation. This method merely assumes Gaussian pulse shapes containing a fixed number of free parameters and the 
Lagrangian form of the partial differential equation is used to obtain the parameter evolution equations. However it is 
a questionable approach because it is based on the a priori assumption that the solution of the PDE has a form which 
remains very close to the chosen ansatz. Accordingly it has to be checked carefully by full numerical simulations of 
the PDE. 

Numerical simulations of the stochastic GP equation with spherically symmetric trap is performed by Crank- 
Nicholson scheme. The absorbing boundary condition is employed to imitate the infinite domain size. This technique 
allows to prevent re-entering of linear waves emitted by the condensate under perturbation into the integration domain. 
We have first checked the variational approach for the unperturbed system. We have done so by inserting the Gaussian 
waveform with the amplitude and width corresponding to a stationary point (as predicted by the variational approach) 
as an initial condition into the PDE JjJJ. We have let the solution evolve in time and we have plotted the results 
in Fig. [SJi. As can be seen the Gaussian ansatz is a good approximation when P is not close to the critical value 
P c - Actually we have found numerically that the critical value for the existence of the BEC is not P c = 0.535, as 
predicted by variational approximation, but P c = 0.459. For P very close to the real value of P c , the Gaussian ansatz 
substantially deviates from the exact solution of the 3D GP equation, as shown in Fig. [U>. 

In a second step we have performed numerical simulations of the GP equation Q driven by a random Gaussian 
white noise r/ with zero-mean and autocorrelation function ~E[rj(t)r)(t')] = a 2 8(t — t'). We do so by choosing randomly 
and independently the value of 77 at each time step. The mean collapse time is calculated as an average over 100 
realizations of random paths along which the width of the condensate evolves from the value corresponding to the 
minimum of the effective potential ao until the value corresponding to its local maximum a± (see Fig. 0. The initial 
wave-form is selected as a Gaussian with parameters predicted by the variational approximation corresponding to the 
stationary state of the condensate. Fig. [7^ represents the collapse time for different values of the parameter P which 
are not too close to the critical value P c . Comparison with the results from numerical simulations of the ODE (J5J 
shows a very good agreement. This demonstrates that the variational approach provides accurate predictions for the 
behavior of the BEC. for small non-linearity, and that the asymptotic analysis carried out in Section llVI holds true 
for the randomly driven GP equation. 

Finally, we have performed numerical simulations of the GP equation JSJl driven by a white noise 77 with a nonlinear 
parameter P very close to the critical value P c = 0.459. For near-critical values of the parameter P the Gaussian 
waveform was found to be not enough accurate. In this case we employed the exact solution of the GP equation 
to initiate random simulations. The exact solution (ground state) of the GP equation is found by imaginary time- 
evolution method as described in [3(j. It is plotted in Fig. The results are plotted in Fig. 03. We can see that 
collapse in the perturbed PDE occurs much earlier than in the ODE model. This shows that the BEC in full GP 
equation is unstable against collapse at near critical nonlinear parameter. A small perturbation can drive the BEC to 
collapse through fluctuations that are not captured by the variational approach. Accordingly, we can state that the 
variational approach provides poor predictions for the behavior of the BEC for critical non-linearity Several reasons 
can explain the departure: 1) the Gaussian ansatz is not correct (see Fig. Et>) - 2) the study of the ODE model shows 
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(a) t (b) r 

FIG. 6: Picture a: Width of the BEC for an initial Gaussian waveform with parameters corresponding to a stationary point 
of the potential U(a). The oscillations are insignificant for small values of P, and become important when P approaches the 
critical value P c — 0.459. At overcritical P the waveform rapidly shrinks (a — * 0), i.e. the BEC undergoes collapse. Picture b: 
Exact solution of the 3D GP equation (solid line) compared with the Gaussian approximation with the same number of atoms 
and P = 0.44. 
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FIG. 7: Mean collapse time calculated from stochastic PDE simulations (solid squares) and compared with the corresponding 
stochastic ODE simulations (open circles). Each mean is computed by averaging over a series of 100 simulations. Picture a: 
Mean collapse time as a function of P for a white noise strength a = 0.3. Picture b: Mean collapse time as a function of a for 
a nonlinear parameter P = 0.44 close to the critical value P c = 0.459. 



that the important parameter in the near-critical case is not the value of P, but the value of the difference between 
P and P c . But the ODE does not capture the correct value of P c , so the error committed in the evaluation of the 
difference P — P c becomes very large when P becomes close to P c . 3) radiation effects become very important, in 
the sense that the waveform is strongly affected, even when the simulations are performed starting from the exact 
numerical waveform plotted in Fig. \^p, so that we feel that it is useless to try to find a more suitable ansatz. In 
this respect, one should add that this result is not surprising because it is well known in nonlinear optics that the 
time-dependent variational approach fails to describe the regime near the collapse 0, • Finally, it is necessary to 
mention that the behavior of the gas close to collapse can be affected by mechanisms that are not included in the GP 
equation, such as inelastic two and three-body collisions |33ll3^|. 



VII. CONCLUSION 



We have considered in this paper a condensate trapped by an external potential generated by a system of laser 
beams in the case of a negative scattering length. We have studied the stability of the metastable BEC against small 
fluctuations of the laser intensity. We have shown that collapse of the BEC occurs whatever the amplitude of the 
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fluctuations after a time which is inversely proportional to the integrated covariance of the autocorrelation function 
of the fluctuations of the laser intensity. The statistical distribution of the collapse time has been computed. The 
dependence of the mean collapse time with respect to the number atoms N has been thoroughly analyzed. We have 
shown that, for N below the critical number of atoms N c , the mean collapse time decreases logarithmically with 
increasing N. As a byproduct of the analysis we have shown that the variational approach is very efficient for the 
analysis of the BEC for a number of atoms N which is not too close to N c , but we have seen that it completely fails 
for N close to N c . 
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